
/*
    Chemistry Formula Balancer

	作者：xkhhdx (xkhhdx [at] gmail [dot] com)
	时间：Aug.10th, 2009

	本程序是一个开源的自由软件，你可以在自由软件基金会的 GNU通用公共许可
	证协议 2.0 的许可下改写或修改本软件。

	修改之后请在你的声明信息中保留本段信息。

	请遵循 GNU GPL v2 协议的条款，具体细节请阅读此协议。

	本文件主要实现了 balancer::balance() 函数。
*/

#include "balancer.h"

#define max( x, y ) ( x > y ? x : y )
#define min( x, y ) ( x > y ? y : x )

#pragma warning (disable:4786)

int solve( double **a, double *b, int n, int *s );

void balancer::balance()
{
	int i = 0, j = 0, n = 0, *s = 0;

	double **a = 0, *b = 0;

	n = matrix_degree;

	b = new double[n];

	a = new double *[n];

	for( i = 0; i < n; i++ )
		a[i] = new double[n];

	s = new int[n];

	for( i = 0; i < n; i++ )
		for( j = 0; j < n; j++ )
			a[i][j] = matrix[i][j];	// 拷贝系数矩阵

	for( i = 0; i < n; i++ )
		b[i] = matrix[i][n];	// 拷贝右边项

#ifdef PROG_DEBUG
	// 输出方程组;
	printf( "\n系数矩阵\n" );
	for( i = 0; i < n; i++ )
	{
		for( j = 0; j < n; j++ )
			printf( "%6.4f\t", a[i][j] );
		printf( " | %6.4f\n", b[i] );
	}
#endif
	// 求解方程组, 并输出方程组的可解信息;
	int ret;

	ret = solve( a, b, n, s );
	if( ret )
		printf( "方程组有解\n" );
	else
		printf( "方程组无唯一解或无解\n" );

#ifdef PROG_DEBUG
	// 输出方程组及其解, 并指出每个变量的唯一性;
	for( i = 0; i < n; i++ )
	{
		for( j = 0; j < n; j++ )
			printf( "%6.4f\t", a[i][j] );
		printf( " | %6.4f\t", b[i] );
		m_coefficient_map[char ( i + 'a' )] = ( int )b[i];

		printf( "\n" );
	}
#endif
	map < char, int >::iterator iter = m_coefficient_map.begin();

	map < char, int >::iterator iter2 = m_coefficient_map.begin();

	cout << "方程解：" << endl;
	for( iter = m_coefficient_map.begin(); iter != m_coefficient_map.end(); iter++ )
		cout << iter->first << " = " << iter->second << endl;

	// 素数表
	int m_prime_table[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };

	// 化简方程，将每个化学计量数除以公因数，得到最简比
	for( i = 0; i <= sizeof( m_prime_table ) / sizeof( int ); i++ )
	{
		int canbedev = 1;

		for( iter = m_coefficient_map.begin(); iter != m_coefficient_map.end(); iter++ )
		{
			if( ( ( iter->second ) % m_prime_table[i] ) != 0 )
			{
				canbedev = 0;
				break;
			}
		}
		if( canbedev )
		{
			for( iter = m_coefficient_map.begin(); iter != m_coefficient_map.end(); iter++ )
				( iter->second ) /= m_prime_table[i];
			// i-- ;
			for( iter2 = m_coefficient_map.begin(); iter2 != m_coefficient_map.end(); iter2++ )
				cout << iter2->first << " = " << iter2->second << endl;

		}
	}


#ifdef PROG_DEBUG
	{
		std::list < chemistry_formula >::iterator iter;
		std::list < chemistry_formula >::iterator iter2;

		int t = 1;

		for( iter = reacts.begin(); iter != reacts.end(); iter++, t++ )
		{
			cout << m_coefficient_map[iter->m_symbol[0]] << iter->m_formula << ( ( t == reacts.size() )? "" : " + " );
		}

		cout << " = ";
		t = 1;

		for( iter2 = outers.begin(); iter2 != outers.end(); iter2++, t++ )
		{
			cout << m_coefficient_map[iter2->m_symbol[0]] << iter2->m_formula << ( ( t == outers.size() )? "" : " + " );
		}
	}
#endif
	
	// Apr.3rd, 2010 这里有内存泄漏的问题…… 已改正
	for( i = 0; i < n; i++ )
		delete [] (a [n]);

	delete[]a;
	delete[]b;
	delete[]s;
}

#define above_zero_min_int 1E-6
#define abs( x ) ( ( x > 0 ) ? x : - x )
#define swap( p, q ) { temp = p; p = q; q = temp; }

/*
 高斯消去法求解方程
 具体的请参见 高中数学 人教B版 必修三 课本的第 38 页或其他资料
*/
int solve( double **a, double *b, int n, int *s )
{
	int i = 0, j = 0, row_pos = 0, col_pos = 0, ik = 0, jk = 0;

	double mik = 0.0f, temp = 0.0f;

	// row_pos 变量标记行循环, col_pos 变量标记列循环

	while( ( row_pos < n ) && ( col_pos < n ) )
	{
		// 选主元
		mik = -1;
		for( i = row_pos; i < n; i++ )
		{
			if( abs( a[i][col_pos] ) > mik )
			{
				mik = abs( a[i][col_pos] );
				ik = i;
			}
		}

		if( mik < above_zero_min_int )
		{
			col_pos++;
			continue;
		}

		// 交换两行
		if( ik != row_pos )
		{
			for( j = col_pos; j < n; j++ )
			{
				swap( a[row_pos][j], a[ik][j] );
				swap( b[row_pos], b[ik] ); // 区域之外？
			}
		}
		// 消元
		b[row_pos] /= a[row_pos][col_pos];
		for( j = n - 1; j >= col_pos; j-- )
			a[row_pos][j] /= a[row_pos][col_pos];
		for( i = 0; i < n; i++ )
		{
			if( i == row_pos )
				continue;
			b[i] -= b[row_pos] * a[i][col_pos];
			for( j = n - 1; j >= col_pos; j-- )
				a[i][j] -= a[row_pos][j] * a[i][col_pos];
		}
		row_pos++, col_pos++;
	}

	// 各个变量解的唯一性
	// 对于每一行a[i] , 若该行是单位向量e[j] , 则解x[j] 唯一
	// x[j] = b[i] , 所以x不另开辟空间
	for( i = 0; i < n; i++ )
		s[i] = 0;
	for( i = 0; i < row_pos; i++ )
	{
		for( jk = 0; jk < n; jk++ )
		{
			if( abs( a[i][jk] - 1 ) < above_zero_min_int )
			{
				s[jk] = 1;
				break;
			}
		}
		for( j = jk + 1; j < n; j++ )
			if( abs( a[i][j] ) > above_zero_min_int )
				s[jk] = 0;
	}

	for( i = row_pos; i < n; i++ )
		if( abs( b[i] ) > above_zero_min_int )
			return 0;
	return 1;
}
